Quantum Monte Carlo simulations of the half-filled two-dimensional Kondo lattice 

model. 
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The 2D half-filled Kondo lattice model with exchange J and nearest neighbor hopping t is con- 
sidered. It is shown that this model belongs to a class of Hamiltonians for which zero-temperature 
auxiliary field Monte Carlo methods may be efficiently applied. We compute the staggered moment, 
spin and quasiparticle gaps on lattice sizes up to 12 x 12. The competition between the RKKY 
interaction and Kondo effect leads to a continuous quantum phase transition between antiferromag- 
netic and spin-gaped insulators. This transition occurs at J/t = 1.45 ± 0.05. 
PACS numbers: 71.27.+a, 71.10.-w, 71.10.Fd 



The Kondo lattice model (KLM) describes a band of 
conduction electrons interacting with local moments via 
an exchange interaction J. This model is relevant for the 
understanding of heavy electron materials Q,^] . The na- 
ture of the ground state results from competing effects. 
The polarization cloud of conduction electrons produced 
by a local moment may be felt by another local moment. 
This provides the mechanism for the Rudcrman-Kittel- 
Kasuya-Yosida (RKKY) interaction || with effective ex- 
change J e ff(q) cx — J 2 Rex(q, w = 0), to) being the 
spin susceptibility of the conduction electrons. On the 
other hand, the same polarization cloud may form a sin- 
glet bound state with the local moment. In the sin- 
gle impurity case, this happens at the Kondo temper- 
ature Tjc oc e^e~ 1 / ,7JV ( e/ ' where et is the Fermi energy 
and N(ef) the density of states Comparing energy 
scales, the RKKY interaction dominates at small J and 
Kondo effect at large J. Thus a quantum transition be- 
tween magnetically ordered and disordered phases is an- 
ticipated. 

The KLM we consider is written as 
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Here i runs over the L 2 -sites of a square lattice, 
corresponds to nearest neighbors, cl creates 
a conduction electron with z-componcnt of spin a 
on site i and periodic boundary conditions are im- 
posed. Si = (1/2) £^,4^'/^, and SS = 

(1/2) £ CT a , cl (Ja.a'C-i a i with a the Pauli matrices. A 
constraint of one fermion per /-site is enforced. At 
J/t <C 1 this model maps onto the periodic Anderson 
model II (PAM) at strong coupling M . Quantum Monte 
Carlo (QMC) methods constitute an efficient tool for 
the study of the PAM in various dimensions I 0. The 
one-dimensional version of the KLM has been extensively 
studied p]] [. In particular, at (n) = 2 (half-band filling 
or one conduction electron per local moment) the Kondo 
effect dominates at all values of J/t. In two dimensions, 
variational Monte Carlo methods [Q as well as series 
expansions around the strong coupling limit |13| support 



the existence of a critical point. The aim of this paper is 
to go beyond the above approximative approaches. We 
show how to efficiently simulate the KLM with the pro- 
jector QMC (PQMC) algorithm This method 
yields exact zero temperature results and is free of the 
notorious sign problem at half-band filling. This stands 
in contrast to previous approaches Jl6| | which generate a 
sign problem even at (n) = 2. Using this algorithm, we 
map out the phase diagram at (n) = 2. 
Our starting point is the Hamiltonian: 



J 



H = -tY cl e, --Y 
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with nl 



fl ff . This Hamiltonian has all the prop- 
erties required for an efficient use of QMC methods. To 
avoid working with continuous fields, we use the approx- 
imate Hubbard Stratonovitch (HS) transformation in- 
troduced in Ref. 0,^8) to decouple the J-term. This 
transformation introduces systematic errors of the order 
(Ar) 3 , where Ar corresponds to an imaginary time step. 
Since this order is higher than the systematic error pro- 
duced by the Trotter decomposition, it is negligible. As 
for the Hubbard term, we have found it essential to use 
Hirsch decomposition in terms of Ising spins which couple 
to the density rather than to the z-component of the mag- 
netization fl9|| . Although this forces us to work with com- 
plex numbers, it conserves SU (2)-symmetry for a given 
HS configuration. As argued in Ref. p0| , this provides an 
efficient algorithm for the calculation of imaginary time 
displaced spin-spin correlation functions pl| from which 
we will determine the spin gap. The ground state of H 
(||), l^o), is obtained by projection. A trial wave func- 
tion, \iPt)i required to be a single Slater determinant 
and non-orthogonal to the ground state is propagated 
along the imaginary time axis till convergence is reached 
|l4],[l5]]. With the above HS transformations and appro- 
priate choice of \^t) [0, particle- hole symmetry leads 
to the absence of sign problem at (n) =2. 
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FIG. 1. Momentum distribution for the /- and c-electrons. 

(a) k = k(l, 1) and (b) k — fc(l,0). Within our precision, 

n/(k) = 1 for all considered values of J/t. To achieve this we 

have chosen values Uf/t ranging from Uf/t = 0.5 (J/t = 0.2) 

to Uf/t = 2 (J/t = 1.6). As J/t grows, n c (k) becomes 

smoother. 

The relation of the above model to the KLM model is 
seen by rewriting Eq. as 

H=-t y ct c? + jy&M- 

jV (d d h „/?„ + h.c. n ) + 

i 

J f r&nl — n5 — nl ) + 

i 

tf/Etet-^Hu- 1 / 2 )- (3) 

i 

with = T"! fJ ff and nS = V* ct o . It is inxpor- 
tant to notice that 

t^X> ~ 4r )(1 ~ n li } + n lAi ] = °- (4) 

i 

The number of doubly occupied and empty /-sites is a 
conserved quantity. Denoting by P n the projection onto 
the subspace with n doubly occupied and empty /-sites 
one obtains: 



HP a = H KLM -N(U f /A + J). 



(5) 



Thus, in principle, it suffices to consider a trial wave 
function satisfying PqI^t) = \^t) to ensure that 
exp(-6ff)|* T ) = exp(~eH KLM )\^T}- The coupled 
constraints: Po\^t) = \^t) and \^t) is a Slater de- 
terminant forces us to choose SI' z \^t) = i^l^r) thus 
breaking SU (2)-spin symmetry. Since the KLM con- 
serves total spin, this symmetry has to be restored by 
the imaginary time propagation. When the energy gap 
to the first excited spin state is small - as is certainly the 
case when long-range magnetic order is present- restor- 
ing this symmetry is extremely expensive. To avoid this 



problem and since the ground state of the KLM at half- 
filling on a bipartite lattice has S = 1 22 2j| , we choose 
a spin singlet trial wave function. During the imagi- 
nary time propagation P n \^T) will be suppressed by a 
factor e ~ nAE ® in comparison to Pq\^t)- In two lim- 
iting cases, we estimate AE ~ Uf/A for J/t <C 1 and 
AE ~ 3Z7//8 for J/t, J/U f > 1. To confirm that we are 
well in the Pq subspace we plot in Fig. p] the single par- 
ticle occupation number nl = (^f \ J2a ft fk o-l^o) w ^ n 

f% a = (l/£) J2j etk J /j cr- Our rc sults are indistinguish- 



able from 



il = 1 which leads to 



a property which may be only realized if Pol^o) = 
|*o)- Owing to Eq. (pi) |\I>o) is nothing but the ground 
state of the KLM. 

We now discuss the phase diagram at (n) = 2 and 
start with the spin degrees of freedom. To estab- 
lish long-range magnetic order, we compute the quan- 
tities S a (r) = |(5 Q (f) • S a {6)) as well as its Fourier 
transform: S a (q) — ^ F e l ^' r S a {r). We consider sep- 
arately the conduction (a = c) and localized (a = 
/) electrons. Long-range antiferromagnetic order is 
present when limj^oo S a (L/2, L/2) = limi^oo S a (Q = 
(?r, 7r))/L 2 takes a finite value. Fig. || plots both above 
quantities versus 1/L for the /-electrons. For lattice sizes 
ranging from L = 6 to L = 12 the QMC data extrapolates 
linearly to a finite value for J/t < 1.45. Similar results 
are plotted in Fig. || for the conduction electrons. The re- 



sulting staggered moment to" = ylimL^oo S a (Q)/L 2 is 
plotted versus J/t in Fig. §. At J/t = 0.2, m{ = 0.557(3) 
- a value much larger than for the Hciscnbcrg model: 
mf = 0.3551(3) p4| . In contrast, m c s is small at small 
values of J/t - m c s = 0.072(6) at J/t = 0.4. In com- 
parison, the half- filled Hubbard model at U/t = 4 leads 
to to s ~ 0.2 0. At the mean-field level, the behav- 
ior of m( ,c at weak coupling may be captured by the 
Ansatz, 01) = (-l) i TO{e 2 /2 (m{ < 1) which leads to 

Hklm = ~*E(5jV i* c l* + ( J7 ^/ 2 ) ZtX- !) ? ^ ' 
Minimizing the free energy with respect to rh( yields: 
to{ = 1. The conduction electrons are thus subject to a 
staggered field of magnitude cx J. Since RexiQ, w = 0) is 
singular, this immediately leads to long-range magnetic 
order with m c s oc ( J/t) ln 2 ( J/t) for J/t < 1. The behav- 
ior of m's'f bears some similarity with the QMC data (see 
Fig. ||). At larger values of J/t the Kondo effect destroys 
magnetic order. Both m{ and m c s scale continuously to 
zero as J/t approaches J c /t ~ 1.45 (see Fig. @). 
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FIG. 2. Spin-spin correlations for the /-electrons versus 
inverse linear length L for several values of J/t. The solid 
lines are least square fits to the form a + b/L. The symbol at 
l/L = corresponds to the extrapolated value. As apparent, 
both S f (L/2, L/2) and S f (Q)/L 2 scale to the same value. 
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FIG. 3. Same as Fig. ^| but for the conduction elec- 
trons. In order to satisfy the relation lim_t_ >00 S* {L/2, L/2) = 
limL_,oo (Q)/L 2 we fit the data to the form a + b/ L + c/ L 2 . 

Once long-range magnetic order is destroyed (J/t > 
1.45) the ground state is expected to evolve smoothly to 
the strong coupling limit, J/t 3> 1. In this limit, \^o) 
is given by a direct product of singlets on the f-c bonds 
of an elementary cell. Starting from this state, a triplet 
excitation has a dispersion relation (up to second order 
in t/J) A^\q) = J - ±§f - fe(q) Q. To compute 
A s (q) numerically we consider S(q,r) = ^(^o\S(q, r) ■ 
S(-q,0)\^ ) where S(q) = S f (q) + S c (q) and S(q,r) = 
e TH S(q)e- TH . For rt > 1, S(q,r) cx exp(-rA s (g)) 
with A s (q) = E (S = l,N,q) - E (S = 0,N). Here 
Eo(S,N,q) denotes the ground state energy with total 
spin S, momentum q and particle number N = L 2 . As 



in the strong coupling limit and for all considered val- 
ues of J/t the spin gap A s = mm.^A s (q) = A S (Q) with 
Q = (7r,7r). Fig. |4|a plots the raw data from which we 
obtain the spin-gap and Fig. ||b, A s versus l/L. A lin- 
ear extrapolation to the thermodynamic limit leads to 
the results plotted in Fig. @. Within our accuracy, the 
value of J/t for which long-range magnetic order van- 
ishes corresponds to the value of J/t where the spin-gap 
vanishes. 

r/t 
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FIG. 4. (a) In S(Q,t) versus rt for J/t = 1.6. The solid 
lines correspond to least square of the tail of S(Q, r) to the 
form a e _rAs . The thus obtained value of the spin gap A s 
is plotted versus l/L in (b). The solid line in (b) are least 
square fits to the form a + b/L. 

Finally, we consider the quasiparticle gap. As appar- 
ent from the single particle occupation number, n c (k), 
(Fig. |l|) the quasiparticle gap grows continuously with 
growing values of J/t. To obtain an accurate estimate of 
this quantity, we compute (^ol Yl<? c j: ( T ) c k J* ) wmcn 
scales as e~ TAqp ^ when rt >■ 1. Here, A qp (k) = 
Eq(N) — Eq(N — 1, k). In the strong coupling limit and 
to first order in t/J, A^J (k) = 3 J/4 — e(k)/2 and thus 
takes a minimum at k = (tt,tt). Our QMC results for 
values of J/t ranging from J/t ~ 0.4 to J/t = 2 are con- 
sistent with: A qp = mmj:A q p(k) = Aqp(-K,ir) |25|] . The 
size scaling of A qp is presented in Fig. [5] and the extrap- 
olated value is plotted in Fig. [^versus J/t. As apparent, 
A qp remains finite and evolves smoothly through the 
quantum transition. In the above discussed mean-field 
approach based on the Ansatz (Si) = (-~l) % m{e z /2, the 
quasiparticle gap scales as J/4 in the small J/t limit. 
Such a behavior, equally seen in one-dimension |]j~l| , p6|| , 
is to a first approximation consistent with our data (see 
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FIG. 5. Same a Fig. ^ but for the quantity 
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FIG. 6. Phase diagram of the KLM at half-band filling and 
T = 0. All plotted quantities are extrapolated to the thermo- 
dynamic limit. At J/t — 0.2 we were not able to distinguish 
Trig from zero. 

To summarize, we have presented an efficient auxil- 
iary field QMC algorithm to simulate zero-temperature 
properties of the KLM. At half-band filling where the 
sign-problem is absent we calculated the staggered mo- 
ment, the spin gap, and the quasiparticle gap on lat- 
tice sizes up to 12 x 12. Our results are summarized in 
Fig. ^|. We observe a continuous quantum phase tran- 
sition between long-range antiferromagnetic and spin- 
gaped phases. This transition occurs at J c /t = 1.45±0.05 
in good agreement with previous approximative results 
fj^ , fl3|| . The quasiparticle gap is finite and evolves con- 
tinuously between both phases. Given that the charge de- 
grees of freedom remain gaped, we expect the observed 



quantum phase transition to belong to the universality 
class of the 0(3) nonlinear sigma model p^,Q. 

A. Muramatsu is thanked for instructive conversations. 
The simulations were carried out on the T3E of the 
HLRS-Stuttgart, as well as on the T90 and T3E of the 
HLRZ-Jiilich. 
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